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I. INTRODUCTION 

For a large class of pattern-forming systems, the essential 
dynamics to be understood and described is that of an inter- 
face between two phases. Mathematically speaking, the prob- 
lem to be solved consists in determining the position of the 
interface as a function of time, i.e., it is a free or moving- 
boundary problem. 

Phase-field models have been established as powerful tools 
for the numerical simulation of this kind of problem. They 
avoid explicit front tracking and are versatile enough to deal 
with topological changes. In a phase-field model, information 
on the interface position is present implicitly, given either as a 
level set of a particular value of the phase field (in two-phase 
models) or by equality of the phase-field values for different 
phases (in multi -phase models), and can be recovered by com- 
putation of the appropriate level set at only those times when 
knowledge of the position is desired. 

A major field of application are solidification problems, 
where diffuse-interface models were developed early on 1 1 , 
2,3] and have seen renewed interest ever since computational 
power increased enough to render their simulation feasible. 
The concept was extended to anisotropic interface properties 
|4J], and first qualitative numerical calculations of dendritic 
growth 0, 01 were followed by theoretical improvement of 
the asymptotics permitting quantitative simulations 001, at 
least for intermediate to large undercoolings. Non-dendritic 
growth morphologies were also simulated, even in three di- 
mensions 1 9]. Generalizations included the description of the 
coexistence of more than two phases 

Additional examples of successful application of the tool 
phase field include the modeling of step flow growth lfl2l [T3] 
and of the elastically induced morphological instability fl4l 
EE, often labeled Grinfeld |Qj] or Asaro-Tiller-Grinfeld 
(ATG) instability ITHl . All of the examples mentioned so far 
dealt with nonconservative interface dynamics, where a par- 
ticle reservoir is provided by either the melt that is in contact 
with the solid or the adatom phase on a vicinal surface. 

Actually, regarding the ATG instability, which is an insta- 



bility with respect to material transport driven by elastic en- 
ergy, interest initially focused on transport by surface diffu- 
sion, which leads to conserved dynamics. This is the case in 
the first article by Asaro and Tiller 1 18], but also in the first 
numerical simulations by sharp-interface continuum models 
fpjTl . preceding computations of the instability under transport 
by melting-crystallization 1 20] . 

The situation reversed, when the phase-field method was 
for the first time employed to compute the ATG instability 
ffl flUl . Here all the early works considered a nonconserved 
phase-field I14lll5lll(il2lll . In fact, there seem to be no publi- 
cations so far treating conserved phase-field dynamics in this 
system. This difference in preferences when modeling either 
on the basis of a sharp-interface model or using a phase field 
may be due to the fact that writing down a nonconservative 
and a conservative model is equally simple in the former case, 
whereas it is less straightforward to write down the conserva- 
tive model within the phase-field approach - correctly - than 
the nonconservative one. The demonstration of this assertion 
will be a major point to be made in the present article. 

This is not to say that phase-field models with a conser- 
vation law for the phase field have not been considered at 
all. Starting from a Cahn-Hilliard equation with a concentra- 
tion dependent mobility, Cahn et al. 1 22] obtained an interface 
equation with the normal velocity proportional to the Lapla- 
cian of the mean curvature. It then looks as if all phase-field 
models with surface diffusion should be derivable on the ba- 
sis of similar considerations. Indeed, comparable models have 
been applied in the simulation of electromigration and voiding 
in thin metal films i23ll24ll . These two models are slightly dif- 
ferent, but the difference is not crucial and all previous models 
seem to suffer from the same problem, to be discussed in the 
following. 

As we shall see, it is quite easy to set up a conservative 
phase-field model. But it is more difficult to obtain the cor- 
rect asymptotics describing surface diffusion as given by the 
desired sharp-interface limit. Past models such as the ones 
presented in 1 22, 23, 24] describe this asymptotics almost cor- 
rectly, but not quite. The purpose of this article is to point out 
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the reasons, to explore an alternative approach and finally, to 
give a model that is asymptotically correct. 

To be definite, I shall consider the case of surface diffusion 
in the ATG instability. The paper is organized as follows. In 
Sec. [n| the sharp-interface model to be approximated by the 
phase-field equations is specified. SectionlHllthen presents the 
standard approach that previously was supposed to reduce to 
the correct limit and pinpoints the oversight in hitherto exist- 
ing asymptotic analyses. An alternative approach is presented 
in Sec. II VI failing for complementary reasons. By an appro- 
priate combination of the ideas from both approaches, I will 
then arrive at a phase-field model in Sec.|V]that may lack el- 
egance but has the virtue of giving the correct asymptotic be- 
haviour and doing so at a lower effective expansion order than 
the preceding models. All derivations are given for the three- 
dimensional case. While this is slightly more involved than 
with the two-dimensional model, it is expected that future in- 
teresting applications will be three-dimensional. Phase-field 
models for conserved dynamics considered so far have been 
restricted to two dimensions. Moreover, the tensorial nature 
of the elasticity part of the model makes a generalization from 
two to three dimensions less straightforward than with scalar 
or vector equations. Section lVTl contains some concluding re- 
marks. A number of mathematical details from differential 
geometry as well as the asymptotics of the elastic part of the 
model, which is unproblematic, are relegated to various ap- 
pendices. 



II. SHARP-INTERFACE MODEL FOR THE ATG 
INSTABILITY UNDER SURFACE DIFFUSION 

One of the more convenient descriptions of the ATG insta- 
bility starts from an expression for the local chemical potential 
difference between the solid and the second phase [liquid, gas 
(vapour), or vacuum] at the interface. In three dimensions, 
this is given by l25ll 
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In this equation, the material parameters are Young's modulus 
E, the Poisson number v, meaning that we assume isotropic 
elastic properties, the surface tension y, assumed isotropic as 
well for simplicity, and the density p s of the solid phase. The 
o-jj are the (cartesian) components of the stress tensor, <x„„ is 
its normal component at the interface, K\ and k 2 are the two 
principal curvatures of the interface at the point considered. 
Their sum is commonly denoted as mean curvature. 

In writing equation Q, we have neglected a gravitational 
contribution to the chemical potential 1 16] that usually is 
unimportant in situations, where surface diffusion governs the 
dynamics. For example, in the self-organized formation of 
patterns under strain built up during molecular beam epitaxy 



of semiconductor heterostructures, the strain energy induced 
by lattice mismatch exceeds any potential energy in the grav- 
itational field by several orders of magnitude. 

To determine the stresses needed in the calculation of the 
chemical potential, one has to solve the elastic equations. 
These are given by the mechanical equilibrium condition of 
vanishing divergence of the stress tensor, 
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where volume forces such as gravity are assumed to be absent 
and we have restricted ourselves to the static limit of elasticity 
(hence there are no sound waves), plus a constitutive equation, 
for which linear elasticity, described by Hooke's law 



■« = rh( Uij + \ ^iy Tj Ukk 5i ) - po6i J ' 



is assumed, in which 



1 / dui duj 
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(4) 



is the (small-)strain tensor. It is convenient to choose the 
hydrostatic strain state assumed by the solid under the equi- 
librium pressure po of the non-solid phase as reference state, 
hence the term -po#y (Po ma Y °f course be zero). This simply 
means that we are measuring displacements of material points 
starting from their position in this state. The issue of refer- 
ence states and their appropriate choice is discussed at length 
in 01. 

The elastic equations must be supplemented with bound- 
ary conditions, of which only the boundary conditions at the 
interface between the two phases need to be specified here: 



0. 



(5) 



that is, the normal component of the stress is equal to the neg- 
ative external pressure, while its shear components along two 
independent tangents ti and t 2 vanish. By assuming a contin- 
uous component <r n „, the pressure jump at a curved interface 
due to capillary overpressure is neglected, an approximation 
that is well-justified under a wide range of conditions. For a 
discussion of the validity of this approximation as well as the 
others mentioned, see 112611 . 

Once the chemical potential has been computed, the stabil- 
ity of the interface can be assessed, but in order to determine 
its evolution, some dynamical law governing its motion must 
be stated. If a particle reservoir is present and the interface 
is rough, it is natural to assume linear nonequilibrium kinet- 
ics. The driving force then is the chemical potential difference 
itself, and the normal velocity v„ of the interface will be pro- 
portional to it: v„ = -k v S/u, where k v is a mobility and the 
normal points from the solid into the second phase. For ma- 
terial transport by surface diffusion, the driving force is the 
gradient of the chemical potential along the surface instead, 
producing a surface current j cc -V s 5/u (V, is the surface gra- 
dient), which leads to a dynamical law of the form 



= D s A s 6p , 



(6) 
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where A s is the Laplace-Beltrami operator on the surface and 
D s a diffusion coefficient, assumed constant here. 

Equations Q through (jfji constitute the continuum model 
for the ATG instability with transport by surface diffusion, 
to which a phase-field model should converge in the limit of 
asymptotically small interface width. 



III. SCALAR-MOBILITY PHASE-FIELD MODEL 

Before discussing the structure of previous phase-field 
models attempting to capture surface diffusion dynamics, let 
us briefly recall the phase-field model for nonconserved order 
parameter 0. This can be written 1 1 611 



at 



= ^{ V V - ^[2/'W + ^h'(<f>)V e i]} (7) 



where /(0) = (p 2 (l - <p) 2 is the usual double-well potential 
describing two-phase equilibrium and h{<p) = 2 (3 - 2(f>) may 
be interpreted as the local solid fraction, varying between 1 in 
the solid phase (i.e., for 0=1) and in the non-solid phase 
(0 = 0). Primes denote derivatives with respect to the argu- 
ment. h((p) has the convenient property that it has extrema at 
= and 0=1, the location of the minima of the double-well 
potential. V e \ is given by 



y el = G y,4 + ^-i(2 
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where G = £/[2(l + v)] is the shear modulus or first Lame 
constant, A = Ev/[(1 + v)(l - 2v)] the second Lame constant, 
and A is the bulk modulus of the non-solid phase (if this is 
vacuum, A = 0). 

In the diffuse-interface formulation, the elastic equations 
take the following form 



dxj 

&ij = h(4>) o- u - [1 - h((p)] p 5 i} , 



(9) 



where cr,y is given by and p - po-A^k u kk- 

The standard approach to a phase-field description of sur- 
face diffusion, as proposed in 1E2II23H24I1 . is then to prepend 
the right hand side of Eq. {7) with a differential operator cor- 
responding to the divergence of a gradient multiplied by a 
phase-field dependent mobility, i.e., Eq. becomes replaced 
with 



-Z = V-MV-^(V 2 0,< 



= -e 2 V 2 + 2f '(</>) + — /i'(0)V e , , 
37 



(10) 



where M is a scalar function of either 1 22, 24] or eV0 |23 
chosen such that the mobility tends to zero far from the in- 
terface: M(0, eV0) — > for — > and — > 1. The elastic 



equations remain unchanged, i.e., while the nonconservative 
model is given by Eqs. @ through the conservative model 
is given by Eqs. (|8j through dlOi . 

At this point, a few remarks are in order. First, the field 
is the density of a conserved quantity by construction, since 
the right hand side of is written as a divergence. This is 
true for any (nonsingular) form of the mobility. Second, 5p 
becomes zero for — > and — > 1, meaning that there is 
no diffusion in the bulk anyway. One might therefore wonder 
whether it is really necessary to choose a mobility that goes 
to zero in the bulk. The conservation law plus the absence of 
diffusion far from the interface should suffice to restrict trans- 
port to diffusion along the interface. In fact, we shall see that 
essentially the same asymptotic results are obtained no matter 
what the form of M, the only conditions to be imposed being 
positivity (for almost all values of or V0) and boundedness. 
It is just easier to derive them if it is in addition assumed that 
M vanishes in the bulk. On the other hand, it will turn out that 
if a restriction imposed by the asymptotics is removed (or not 
yet satisfied in the temporal evolution of the system), M has 
to decay sufficiently fast inside the bulk for the limit to make 
sense. This may be relevant for the behaviour of the model 
before it reaches its asymptotic state. 

Finally, the issue at present is not so much whether the dy- 
namics is conservative but whether it does reduce to the sharp- 
interface model of Sec. [H] in the limit of an asymptotically 
vanishing interface thickness. To investigate this, we have to 
explicitly carry out the asymptotic analysis. 



A. Local coordinate system 

The basic idea of the analysis is to expand all dynamical 
quantitities in terms of the small parameter e describing the 
interface thickness, to solve for the phase field and to use the 
solution to eliminate its explicit appearance from the equa- 
tions. To this end, the domain of definition of the field is di- 
vided into an outer region, where gradients of the fields can be 
considered to be of order one and an inner region (close to the 
interface), where these gradients are of order 1/e. The expan- 
sion in powers of e is rather straightforward in the outer do- 
main, Eqs. (JSJi through fllOi can be taken directly as a starting 
point. As to the inner domain (and its matching with the outer 
region), it is useful to first transform to coordinates adapted 
to its geometry. Therefore, a coordinate r is introduced that 
is orthogonal to the interface (defined as the level set corre- 
sponding to (p(x,y,z,t) = 1/2) and two coordinates s and q 
that are tangential to it. r is the signed distance from the inter- 
face and will be rescaled by a stretching transformation r - ep 
to make explicit the e dependence for the expansion. Inner and 
outer solutions must satisfy certain matching conditions due to 
the requirement that they agree in the combined limit e — > 0, 
p — > +00, r — > 0. These conditions are given in Add. 1X1 

To obtain a basis for our local coordinate system, we first 
write 

r = R(s,q) + r w(s,q) (11) 
where r is the position vector of a point near the interface, 
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R the position of the interface itself, and n the normal vector 
on it (oriented the same way as in the sharp-interface model, 
i.e., pointing out of the solid). Next, we require the coordi- 
nates s and q parametrizing the interface to be both orthog- 
onal and arclength-like, meaning that the tangential vectors 
ti = dR/ds, t 2 = dR/dq are orthogonal and that they are unit 
vectors. It is shown in the App.|B]that the requirement of both 
coordinates being arclengths is less innocent than it seems - 
it implies that the coordinate lines are geodesies, which is a 
strong restriction. (On a planar surface, only cartesian coordi- 
nate systems satisfy this condition.) It may be worthwhile to 
point out a difference between the two- and three-dimensional 
cases here. In two dimensions, the local coordinate system 
can essentially always be chosen global in the tangential co- 
ordinate, for which we may take the arclength. Locality is 
only necessary in r (this coordinate may turn singular for a 
curved interface, when the distance from the latter reaches the 
radius of curvature). In three dimensions, an interface with 
a nontrivial topology does not permit a description by a sin- 
gle coordinate patch, so in general the coordinates s and q 
have to be local, too. However, for the asymptotic analysis, 
we never need to compare far-away points along the interface, 
so a coordinate patch that is valid in the neighbourhood of a 
point R(s(), qo) is quite sufficient, and we can choose, without 
loss of generality, a geodesic coordinate system, despite the 
fact that this will not necessarily be global even on a curved 
interface with the topology of the plane. 

Given the coordinates, it is a trivial matter to write down a 
coordinate basis 



£,.= — = n(s, q) , 
or 

g = * - ^ +r ?H 
s ds ds ds 
= dr _ dR +r dn 

q dq dq dq 



(1 +rK X )t x +r Tl t 2 . 
■ (1 + rK 2 )t 2 - rr 2 ti 



(12) 



The second equality in each of these equations, i.e., the rep- 
resentation of the coordinate basis in terms of the (right- 
handed) orthonormal basis (n,ti,t2) is derived in Apd.IbI k\ 
and k 2 are the normal curvatures of the curves q = const, and 
s = const., respectively. Because the geodesic curvatures of 
these are zero, they are identical to the total spatial curvatures. 
Their sign is chosen such that a protrusion has positive curva- 
ture, i.e., a sphere would have positive curvatures, opposite 
to the traditional definition in mathematics. t\ and t 2 are 
the torsions of these curves, and in Add. 151 it is shown that 
ri = -t 2 = t. 

Note that the basis vectors S s and & q are guaranteed to be 
orthogonal to each other only for r — 0, i.e., on the interface. 
It would certainly be possible to define coordinates that are or- 
thogonal in a whole neighbourhood of the interface point un- 
der consideration, but the additional effort would not be justi- 
fied by the minor simplifications afforded. What will become 
important in the following, however, is that £, is orthogonal to 
the two other basis vectors. Moreover, has unit length for 
all values of r, while the lengths of the two other basis vectors 
vary off the interface (i.e., for r + 0). 

For many of the calculations to follow it is convenient to 



stick to the nonnormalized coordinate basis most of the time 
and to use the orthonormal basis only at the very end, when 
the dependence on simply interpretable quantities such as the 
curvatures is to be exhibited. This enables rather compact rep- 
resentations of most expressions. From d!2i . we first obtain 
the metric coefficients g a p = & a &p, where a,fi e {r, s, q). The 
metric tensor then reads 



(gap) = g 

10 ^ 

(1 + rKi) 2 + r 2 T 2 2rr + r 2 {K\ + k 2 )t 

2rr + r 2 {K\ + k 2 )t (1 + rK 2 ) 2 + r 2 r 2 

its determinant is 

g = detg = [(1 + + rK 2 ) - r 2 T 2 f , 



(13) 



(14) 



and the contravariant components of the metric tensor are ob- 
tained as 



1 




^ 

(1 + rK 2 ) 2 + r 2 T 2 -2rr - r 2 (Ki + k 2 )t 
■2rr - r 2 {K\ + k 2 )t (1 + rK\) 2 + r 2 j 2 



(15) 



From now on, we use the Einstein summation convention 
for pairs of covariant and contravariant indices. The vectors 
of the reciprocal basis are obtained from & a = g^&p. 



S r = Vr = n(i, q) , 

8 s = Vs = ((1 + wc 2 )ti - rrt 2 ) , 

fi ? = V« = ^ = ((l+r*i)1a-rrti) . 

The gradient and divergence read 
V = 6 a d a , 

V-A=-L,9 ff (v?^%) • 



(16) 



(17) 
(18) 



In the following, we will denote inner quantities by the up- 
percase letter corresponding to the lowercase letter denoting 
the outer quantity, whenever this is meaningful. 

Since the interface will move in general and the coordinates 
r, s, q are defined with respect to the interface, there is also a 
transformation rule for the time derivative: 



d,f(x, y, z, t) = d,F(r, s, q, t) - vVF(r, s, q, t) . 



(19) 



where v(s, q, f) is the interface velocity. Equation dl9l ex- 
hibits that the time derivative in the comoving frame is a mate- 
rial derivative. In order to formulate the matching conditions 
concisely, we will occasionally also write the outer fields as 
functions of the variables r, s, and q (without changing their 
naming letter, thus in this case adhering to the physicists' con- 
vention of using a letter for a physical quantity rather than a 
mathematical function). 
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B. Inner equations 

To render the scales of the different terms more visible, we 
rewrite Eqs. d!7t and dl 81 in the (still exact) form 



V= - nd p +8 a d a , 



(20) 



V • A = -L |I dp V?A r + «9 S V&*%) , (21) 

where sub- and superscripts carrying an overbar are taken 
from the set {s, q] only in the summation. 

Assuming, without loss of generality, that the tangential ve- 
locity of the interface vanishes, Eq. fllOi takes the following 
form 

1 1 , 

<9 r O - -v„<9 p <D = V-MV^^(V 2 0),<D), 



6 P ~ ' ""e 2 



^(V 2 0),<D) = - dp ^d p <S> 



ys 

+2fm + ^h'(<I>)V e i , (22) 
3y 



with 



V ■ MV = 1 4= fl P <8 Md p + 4= fl « <88 a ~ p Mdp . (23) 
e 2 V# V? 

Hence, the leading term of the inner equation (12 2 1 with the 
differential operator given by J23l > is of order e~ 4 . 
The mechanical equilibrium condition becomes 



(24) 



C. Expansions, matched asymptotic analysis 

To solve the outer and inner equations successively, we ex- 
pand the fields in both the outer and inner domains in powers 
of e 

</>(x,y,z,f) = </> (0 \x,y,z,t) + e^ m {x,y,z,f) 

+ e 1 <p (2 \x,y,z,t)..., (25) 



iijj{x, y, z, t) = uf){x, y, z, t) + euf){x, y, z, t) 
+ e 2 u (2) (x,y,z,t)... , 



(26) 



and 



3>(r, s, q, t) = <b {0) (r,s,q,t) + e<& m (r,s,q,t) 

+ £ 2 <& (2) (r,s,q,t)..., (27) 



U ofi (r, s, q, t) = C/g(r, s, q, i) + eU { "{r, s, q, t) 
+ £ 2 U^(r,s,q,t).... 



The transformation from cartesian tensor components (sub- 
scripts to curvilinear ones (sub- and superscripts a,/3) is 
carried out using standard rules of tensor calculus. Details 
can be found in App. From the expansions of the strain 
fields, similar expansions (not written out here) follow for the 
stress fields. Moreover, the basic field variables for which the 
elastic equations, the Lame equations, ultimately constitute a 
closed set, are the displacement vectors m, or U a . A pecu- 
liarity of phase-field formulations of elasticity is then that the 
fields are independent of p fl^l . i.e., continuous across 
the interface, because components of the strain tensor con- 
tain derivatives d p U a /e, but their leading order is e°, not e~ l . 
Again, details are presented in App.lCl 

The asymptotic analysis of the model may be divided into 
two natural subproblems: analysis of the equation for the 
phase field and analysis of the mechanical equations; the latter 
depend on the former only via the solution for the phase field. 
Since the focus of this article is on how to capture surface dif- 
fusion dynamics correctly, only the first subproblem will be 
treated in the main text, a strategy that will hopefully facilitate 
the application of the results to other models, such as phase- 
field simulations of electromigration, for which they should 
possess relevance, too. The derivation of the sharp-interface 
limit of the mechanical equations is a straightforward exten- 
sion of the two-dimensional case, albeit with one or two addi- 
tional subtle points. Because the three-dimensional case does 
not seem to have been discussed in the literature before, this 
is done in App. [C] Another motivation to carry this part of 
the analysis out explicitly is to present, as a case in point, an 
operational phase-field model for a complete physical system. 



1. Leading order 



The leading-order outer equation for (f> is 
V ■ MV/(0 (O) ) = , 



(29) 



(28) 



which is to be supplemented with the boundary conditions 
= 1 and (O) = at infinity in the regions where the sys- 
tem is solid and non-solid, respectively. If we regard ( I29t as 
a partial differential equation for the function /'(^ <0) ) (rather 
than for (O) itself), this boundary condition translates into 
/'(0 (O) ) — > as |r| — > oo, which may be seen immediately 
from the explicit form of f'{<p), given in App. [D] The new 
boundary condition is valid everywhere at infinity except pos- 
sibly in a region with size of order e. For general M{<p, eV</>), 
the partial differential equation J29l > is of course nonlinear. 
Nevertheless, it can be shown to have the unique solution 
/'(0 (O) ) = 0, if M is positive everywhere, except possibly on a 
set of measure zero. 

To see this, multiply Eq. \29\ by /'(0 (O) ), integrate over all 
of space and use Gauss's theorem: 

= J^ <i 3 x/'(0 (O) )V • MV/'(0 (O) ) 
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where the 0(e) stands for the surface integral at infinity. If M 
is positive almost everywhere, we immediately get f'((f> (0 ^) = 
const., and the boundary conditions require the constant to 
be zero. This conclusion remains of course unchanged, if M 
becomes zero only when (O) is zero or one - a standard choice 
His M(<p)cccf> 2 (l-(f>) 2 . 

Hence, the unique solution to the leading-order outer prob- 
lem is, if we now consider it an equation for </> (0) again, </> (0) = 
1 in O and (O) = in Q + , where Q T are those regions of 
space, separated by the interface(s), in which limir^oo (O) = 1 
and 0, respectively. The solution (O) = i, still possible for the 
equation interpreted as an equation for /'(0 <O) ), is excluded by 
the boundary conditions for <p^°\ (This argument presupposes 
that we have no domains that are not connected with infinity. 
For the interior of a closed interface, the solution (O) = 1/2 
would have to be excluded by a stability argument or by mak- 
ing reference to initial conditions.) 

It is then seen by inspection that the outer equation is indeed 
solved to all orders by the solution under discussion (because 
h'(0) = h'(l) = 0). Therefore, we have = 0, <f> (2) = 0, 
which provides us with partial boundary conditions for the 
inner solutions <D (1) , 4> <2) , and so on (see App. lA"t. More- 
over, only the inner problem needs to be considered beyond 
the leading order. 

Because g — 1 + 0(e), the leading-order inner problem be- 
comes [see Eqs. J22l > and J23H 

d p M(^)d p [d pp O (0) - 2/'(O (0) )] = , (31) 

which can be integrated once to yield 



d p [d pp ^-2f'(^)] = 



ci(s, q) 
M(0>(°>) ' 



(32) 



where c\(s, q) is a function of integration. It is here that 
we have to follow different lines of arguments, depending on 
whether M approaches zero for p — > +oo, which is the case 
for the mobilities assumed in 1 22, 12411 . or whether it is just 
a bounded (and possibly constant) function of (p. In the first 
case, we may immediately conclude c\ = 0, because the right 
hand side of J32t must not diverge. In the second case, we ob- 
tain the same result by integrating d32l > first and invoking the 
boundary conditions: 

d pp O (0) - 2/'(<D (0) ) = ci(i, q) f \-dp + c 2 (s, q) . (33) 

Jo M 

Since M is bounded from above and positive, the integral will 

CP 

be larger in magnitude than J l/(sup p M)a!p = p/sup p M, so 
the two terms on the right hand side of Eq. d33> are linearly 
independent. The left hand side approaches zero for p — > ±oo 
[the argument will be made more rigorous below in the dis- 
cussion of <I> (1) ], so both c\ and c 2 must be equal to zero. To 
argue that c 2 is zero in the case where M — > for p — > +oo, we 
can proceed the same way, except that we have already gotten 
rid of the term containing c\, so the right hand side of (1331 is 
c 2 only. 

To summarize, the leading-order inner equation results in 



and this provides us with the solution (D (0) = l /2 ( 1 - tanh p) as 
is shown in Adp.IdI 



2. Next-to-leading order 

The next-to-leading order in Eq. (I22i is the order e~ 3 . Since 
the differential operator in front of the chemical potential is 
of order e~ 2 and the chemical potential multiplied by another 
factor e~ 2 , we must expand dp up to order e. Equation ( 1341 
already tells us that 6f/ 0) = 0, so we obtain 



d p M(¥ 0) )d p Sp (l) = 



from which we get 



d p 6/i 



d\(s,q) 
M(<D(°>) ' 



(35) 



(36) 



As before, we can immediately conclude from this that d\ - 0, 
if we assume M(O (0) ) -> for <£ (0) -> 0, 1. For arbitrary but 
bounded M, we invoke the matching conditions (see Apd.IATi 



lim d p 5p ( ' = d r 5p 

p— >±oo 



(1) - * -4% 







(37) 



to obtain the same result (where for once we have denoted an 
outer quantity by a subscript "out") 

Integrating once more with respect to p and writing out 
6p( l \ we have 

SiP = -3 PP <D (1) - (*i + K2)3 P 0> (0) 

+ 2/"(<D (0) )(D (1) + — h'(^)Vf = d 2 (s,q) . 
3y 



(38) 



Up to this point, there is agreement between this and preced- 
ing asymptotic analyses 12311 . if not in all details of the pro- 
cedure, so at least in the results (assuming the appropriate re- 
placement of the interaction with the electric field in the elec- 
tromigration system by that with elastic strain and simplifying 
our result to the two-dimensional case). 

Let us now try to determine the function of integration 
d 2 (s, q). A priori, there is no reason to use a procedure dif- 
ferent from what we have done before. We know the limiting 
values for p — > ±oo for three of the five terms on the right 
hand side of J38i : lim p ^ ±0O (9 P <D <0) = [which follows from 
either the matching conditions or by inspection of the solu- 
tion lIDSl l. lim p ^ ±0O /i'(<D (0) ) = [because h'(0) = h'(l) = 0], 
lim p ^ ±M d 2 (s, q) = d 2 (s, q) (because d 2 is independent of p). 
Moreover, from the matching conditions, we obtain the limit 
for (D (1) 



P0' (O) (+O) + U '(±O) 



.(!)/ 



<f> (1) (+0) = 

(p -> ±co) . (39) 



3 PP <D (0) - 2/'(O (0) ) = , 



(34) 



The second equality follows from the fact that (O) = or 
(O) = 1, hence its derivative with respect to r vanishes on 
both sides of the interface; the third equality is a consequence 
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of the fact that </> (()) solves the outer equation to all orders and 
hence <p^ = 0. 

With four of the five terms in J38I having a definite limit, 
we may conclude that the fifth must have a limit as well and 
obtain 



lim -3 pp O (1) = d 2 (s,q). 

p— >±co 



(40) 



But if this limit exists, it cannot be different from zero: trans- 
forming to ^ = 1/p, we see that dpp<D (1) = (£ 2 3 f ) 2 <I> (1) , which 

implies the asymptotic behaviour <6 (1) d 2 /2£ 2 — > 0) and 

hence the divergence of <D (1) as -d 2 p 2 /2, \fd 2 * 0. (The same 
kind of argument can be used to show that the left hand side 
of Eq. (1331 goes to zero, even though the matching conditions 
do not provide a direct expression for lim p ^ ±ca <9 P p<E> (0) .) 



The upshot of these detailed considerations is that 
d 2 (s, q) = => 6p m = 0. 



(41) 



Previous treatments of the problem did not enter into these 
considerations. Instead, one of the two following equivalent 
approaches was chosen. Either, Eq. (I38i was interpreted as a 
linear inhomogeneous differential equation for <£ (1) and Fred- 
holm's alternative invoked. Since the appearing linear opera- 
tor 



L = 8 PP - 2/"(<D <0 >) 



(42) 



is hermitean, we know (from taking the derivative of Eq. (I34> 
w.r.t. p) that <9 P $ (0) is a solution to the adjoint homogeneous 
equation. The inhomogeneity of the differential equation must 
be orthogonal to this solution. Or else, Fredholm's alternative 
was not mentioned, the equation was simply multiplied by 
3p<l> (0) , integrated, and it was shown via integration by parts 
that the terms containing disappear. Of course, this is the 
same thing. 

To exploit the approach optimally, we use the relationship 
(9pO (0) = -/!'(<D <0) )/3 shown in App.|D] This way we obtain 
from (ED 



oc 

-r 



((» 



Kl+K 2 + -V. 

y 



{d p &°fdp 



d p <& (0) d 2 (s,q)dp = -d 2 (s,q) 



</j(.s.</) = i(*i+«fc + v2°/y) 



(43) 



(44) 



>(0) 



where we have used jDllH and introduced V ei 

Both Eqs. fl41i and 144> were derived by valid methods, 
therefore they should both hold true. Nevertheless, as we shall 
see shortly, Eq. J41i is a quite undesirable result. This may be 
the deeper reason why it was so far overlooked and only the 
analog of Eq. (1441 derived. When Eq. fl41i is inserted in (1441 it 
imposes a relationship between the elastic state of the material 
(represented by V^) and the curvature. (In models, where the 
interaction term is quadratic in e 12311 - it would even impose 
the restriction k\ + k 2 — 0.) 



3. Higher orders 

To see that the model would indeed work if we did not have 
the restriction (14 1 i . let us consider the equations at the next 
two orders, ignoring for the time being the result d 2 = 0. Since 
both (5p (0) (= 0) and 6p (V) are independent of p, the first term 
of the operator d23t does not produce any contribution from 
these terms in (I22> . and the order e~ 2 equation reads 



d p Md p 6p a) + d 5 g^Mdp6p {0) = . 



(45) 



where we can immediately drop the second term, because of 
£p (0) = 0. After two integrations this becomes 



6n (Z) = e^s,q) 



Jo M 



dp + e 2 (s,q) 



(46) 



If M — » for p — > +oo, we immediately find e\(s, q) = 0. In 
the general case, we use the matching conditions [see (IA6H 



Sp 



(2) 



1 

2 

+ <^outl<-=±0 ■ 



^p 2 d rr 6fxf^\ r=±0 +pd r 8p m 



I outl'-=±0 



(47) 



From Eq. (1101 . we gather that an expansion of 6p ont in pow- 
ers of e will contain three types of terms, the first of which 
have the form V 2 </> (<;) (k = 0,1,...), while the second con- 
tain factors (f> (k) (k = 1,2,...), coming from an expansion of 
f'(<p) or h'((p) about cf> (0 \ and the third include either /' (0 (O) ) 
or h'(0 (O) )vf (k = 0, 1,...)- All of these terms vanish, be- 
cause W = for k > and because /'(0 <O) ) = h'((/) (0) ) = 0. 
This is simply a consequence of the fact that the outer equa- 
tion is solved exactly by (O) = and </> (0) = 1 . The "chemical 
potential" appearing in the phase-field equations needs to be 
related to the true, i.e., sharp-interface chemical potential only 
at the interface. In the outer domain, it is zero. We can then 
conclude from (1471 that e\{s,q) = (of course e 2 (s,q) = 0, 
too, but we shall not make use of this result). 

Given that SpP^ is independent of p, the inner equation at 
order e~' takes the form 

-v„d p ¥ 0) = d p Md p 8p 0) + dag^Md-p dp {l) . (48) 

Because the surface metric g°P = S^a to lowest order in e, we 
find after integration over p (v„ does not depend on p) 

m((D (0 Vp 

00 

+ M(<D (0) )a p ^ (3) |!° co . (49) 

Here we can drop the second term on the right hand side, 
if ]im p ^ ±m M(O <0) ) = 0. Formally setting J ^ M(O (0) )Jp = 
3M*y and using (1441 . we arrive at 

v n = M* (d ss + 8 qq ) [y( M + k 2 ) + Vf] . (50) 

It is shown in App. IClthat Vy is the correct sharp-interface 
limit of the elastic part of the chemical potential in Eq. Q, 



8 



hence ( 15 01 reproduces the sharp-interface limit 0, as the sur- 
face Laplacian in geodesic orthogonal coordinates is given by 

A s =d ss +d qq . (51) 

Finally, M* would be infinite for positive functions M(O (0) ) 
that do not reduce to zero for p — > ±oo; therefore, in the end 
we would indeed have to require M(O <0) ) to decay far from 
the interface, if 6/j (1} were different from zero. In reality, we 
do not just have (15 0> . the equation we want, but in addition 
Eq. (I4H . requiring <5ju (1) = and, consequently 

v„ = . (52) 

Thus we are led to the conclusion that the asymptotics of the 
phase-field model produces the correct sharp-interface limit 
only at equilibrium, i.e., when the normal velocity vanishes. 
Of course, one may ask what will happen, if we prepare the 
system in an initial state, where Eq. J10I requires the phase 
field to change. Then v„ will necessarily be different from 
zero. The answer is that in such a case the phase field will not 
follow its asymptotic dynamics yet. A similar phenomenon 
happens when a phase-field simulation is started with an ini- 
tial interface perturbed by white noise. Since the asymptotics 
of the phase-field equations require curvatures to be smaller 
than 1 /e, the initial stage of the dynamics where larger curva- 
tures are present, will not be governed by these asymptotics. 
But the asymptotic behaviour is sufficiently robust to keep that 
initial stage short. It is tempting to speculate that when con- 
ditions are such that (I52> does not hold yet, the phase-field 
model discussed here will satisfy all less restrictive conditions 
such as ( I50> already. Then the model would be applicable dur- 
ing the period where the influence of condition (14 1 i leading to 
(I52> is still small. However, it should be clear that without 
a theoretical estimate of the error in this not fully asymptotic 
state, the model can hardly be considered quantitative. Condi- 
tion J41i should be expected to have a stabilizing influence on 
the interface, rendering growth of instabilities less rapid than 
in the sharp-interface model. 

At first sight, the deep quench case considered in l22ll 
seems to escape an analogical consequence, because there the 
matching procedure is performed on a finite p domain and 
hence an argument leading to the restriction J4 1 ft appears to 
be missing. Indeed, instead of /(</>) the model has a quadratic 
potential /(</>) = (1 - </> 2 )/2, leading to sinusoidal behaviour 
of <D (0) , which takes its limiting values -1 and 1 for finite p 
(= +n/2) already. However, closer inspection shows that the 
problem has not disappeared, it just takes a hidden form. The 
phase field must be specified beyond the domain of thickness 
r = 0(e) where f(<p) + 0, which means that the definitions of 
/'(<D (0) ) (= -<D (0) ) and /"(<D (0) ) (= -1) have to be extended 
to include the values <t> (0) — ±1 into their domain. To keep 
Eq. (3.10b) of 1 22] valid for arbitrarily large p, it is necessary 
to set /'(±1) = 0, i.e., to impose jump discontinuities of size 
1 at the ends of the cD (0) interval [-1,1]. This is evident, as 
the stationary values of the phase field should correspond to 
minima of the potential. However, Eq. (3.13b) of 112 211 then 
becomes 

p(s, t) = -0$ - *<D< 0) + /"(0> (0) ) <D (1) , (53) 



where the behaviour of the second derivative of the potential 
is essentially given by /"(<D (0) ) = -1 +5((D (0) - 1)-k5(O (0) + 1). 
In any case, with /'(O <0) ) exhibiting only jump discontinuies, 
y"(<j)(0)) canno t have stronger singularities than 5 functions. 
The last term in Eq. d53i then becomes zero for \p\ > n/2, 
because approaches zero continuously forp — > ±n/2 (and 
is identically for \p\ > n/2). Then the total right hand side 
of the equation will approach zero for p — > ±00, again forcing 
the chemical potential p(s, t) to be equal to zero. 

It is instructive to note why the nonconservative model ob- 
tained when dim is replaced with Q does not suffer from a 
similar difficulty. In that model, the velocity is already deter- 
mined at the next-to leading order. Instead of ( I38> . we get 

Ps 



- 2/"(<D (0) ) <D (1) - 




Again we may conclude that all the terms on the right hand 
side go to zero as p is sent to +00. However, this does not lead 
to any constraints, since the left hand side is p dependent now 
and goes to zero as well, satisfying the limit automatically, 
whereas in the surface-diffusion case, it was a function of s 
and q only (^2) that could be concluded to be equal to zero. 
So consideration of the limit does not produce anything new 
here, and the only procedure available to extract information 
on v„ is to use Fredholm's alternative which gives the correct 
sharp-interface limit. 

In the case of the nonconservative model, the introduced 
chemical potential functional is zero in the bulk just as in the 
conservative case, but there are no restrictions on its varia- 
tion near the interface, where it acquires a form tending to 
a 6 function in the sharp-interface limit. In the conservative 
model, this is excluded by restrictions on the derivative of the 
chemical potential with respect to p, meaning that the latter 
must be smooth across the interface. Since it is zero off the 
interface, it is zero on it as well. Due to this reason, the phase- 
field model strictly speaking applies only to the equilibrium 
limit. Far-from equilibrium dynamics is not likely to be cap- 
tured faithfully. 

IV. TENSORIAL MOBILITY 

That the phase-field model given by Eqs. (|8} through JlOl t 
does not quite yield the correct asymptotics may be traced 
back to the fact that the differential operator V-MV, prepended 
to the chemical potential in fllOi . does not reduce to the surface 
Laplacian A s in the asymptotic limit. In fact, the second term 
on the right hand side of Eq. i23\ is the Laplace-Beltrami op- 
erator on the surface, but the first term, containing derivatives 
with respect to p is orders of magnitude larger, being preceded 
by a factor of 1/e 2 . As a consequence, the asymptotics must 
be secured by the full solution of the equation rather than by 
both the operator and the chemical potential converging to the 
desired sharp-interface limits. 
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Realizing this property of the model, it seems natural to 
modify the differential operator via introduction of an essen- 
tially tensorial mobility. Let us denote by 



n = 

m 

the normal on the surface <p = const, (for . 
n = n), then we expect the operator 

V-(l-n:n)V 



(55) 

1/2, we have 
(56) 



to reduce to the surface Laplacian asymptotically. A colon is 
used to designate a dyadic product, so 1 - n : n is a projection 
operator projecting onto the tangential plane of a level set of 
(f>. Introducing the shorthand <f>, a = d a <p, we have V0 = S a (p, a 
and 

V-(l - n : n)V = 



(V0) 2 



^ v g<^, a cU. (57, 



The third expression is obtained from the second applying the 
divergence operator Jl 8I > and renaming a — a, n /3 
in the three pairs of "mute" indices. 

To expand this operator in powers of e in the inner domain, 
we use the fact that g pa = and that 

(VO) 2 = 4"(V3>) 2 , 

e l 

(VO) 2 = e 2 (D,„^(l>, /3 =(D, 2 +e 2 0, s /%, j 8 . (58) 

Inserting this into (I57> and carrying the expansion to formal 
order e , we find first that the order e~ 2 terms (containing two 
derivatives with respect to p) cancel each other. The remainder 
reads 



V-(l -n : n)V = 
1 



■■dpyfg 



(VO) 2 



1 -/ 0, 0,n 



(V(D) 2 



0(e). 



(59) 



and this expression still contains derivatives with respect to 
p. However, ;/ the leading-order solution O (0) depends on p 
only, as it did in the last section, then all the derivatives of 
O with respect to s and q (i.e., with respect to the variables 
marked by an overbar) are 0(e) at least, and Eq. d59t reduces 
to V-(l-n : n)V = l /tfd? yfgg vft df,+0(e), i.e., at leading order 
the operator indeed becomes the Laplace-Beltrami operator 
on the surface. 

This then suggests to replace the phase-field equation ( fTol 
with 



8<h 

— = MV ■ (1 

dt 



n: n)V-U/z(V 2 



(60) 



where 5fi is unchanged from (II Oi but M is a constant mobility 
now. 

In this model, the equation for the velocity would appear at 
the next-to leading order already and take the form 

v„S p O (0) = M (d ss + d qq ) {d pp & l) + fa + K 2 )d p ® (0) 



2/"(<D (0) )(I) (I) -— /7'((D (0) )y^j 
3y J 



(61) 



Because the operators X [defined in Eq. J42H and d ss + d qq 
commute, Fredholm's alternative could be applied as in the 
nonconservative case. This eliminates <t> (1) from the equation 
and produces the correct sharp-interface limit. 

In spite of this enjoyable state of affairs, model J60i fails 
much more miserably than dlOi . The reason is that the zeroth- 
order solution is not unique. In fact, the leading-order outer 
equation 



V-(l-n«»:fi (0) )V/V 0) ) = 



(62) 



is solved by any (differentiable) function (p^ satisfying 
the boundary conditions: obviously we have V/'(0 (O) ) = 
f"(<p (0) W<p {0 \ whence 



n 



(0) . ft(0)r 



n l0| V/'(/») = -n (U, /"(0 (U) )|V^ 



(0) f"fA<»\ 



a(0)i 



= /"W ( W=V/'(f), (63) 



which implies (1 - n (0) : n (()) )V/'(^ (0) ) = for all functions /' 
of^°». 

It can be said that this model fails for reasons complemen- 
tary to those of the scalar model. Whereas we had one equa- 
tion too many in that case, adding a constraint to the desired 
sharp-interface dynamics, now we have one equation too few, 
as there is nothing in the model fixing (O) . If we had the right 
(O) , the tensorial model would work perfectly. In particular, 
it would give the limiting equations at lower order than the 
scalar model, which may have positive implications for the 
numerical cost at a given desired level of accuracy, as will be 
discussed in the conclusions. 



V. MODEL WITH CORRECT ASYMPTOTICS - 
MODIFIED TENSORIAL MOBILITY 

In order to obtain a model not plagued by either of the dis- 
advantages of the two cases discussed, it appears that we have 
to combine ideas from both. While it is certainly desirable to 
have a differential operator that itself approaches the surface 
Laplacian, it should do so only for phase field functions that 
have the correct leading-order profile. 

This can be achieved by modifying 1 — n : n into 



2=1 



, V0 : V(p 



1 



e 2 (V0) 2 



-n : n . 



(64) 



4f(4>) 4/(0) 

If we replace the projection operator in Eq. ( I60> by Q, then 
the outer equation at leading order will have the same differ- 
ential operator as the scalar model with constant M. 
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On the other hand, in the inner domain, we have, pro- 
vided O (0) solves the differential equation jD6> . e 2 (V<E>) 2 = 
O, 2 +0(e 2 ) = 4/(0) + 0(e) [this follows from Eqs. l|58}, IDTl . 
and JD1H . whence Q « 1 — n : n. 

This approximation is accurate up to (9(e) only, which is 
not sufficient, because the order e correction would enter as 
a bothersome additive term in the next-to leading order inner 
equation. 

A better inner approximation to 1 - n : n than just Q is 
provided by a minor modification. Obviously, we have Q = 
1 - ft : ft + (9(e) ft : ft in the inner region. Taking this to 
some integer power m we get, because 1 - ft : ft and ft : ft are 
orthogonal projectors: 



Q m = 1 - ft : ft + 0(^) ft : ft . 
These considerations lead us to make the ansatz 

at e z 
6p(V 2 cf>,<f>) ee -e 2 V 2 + 2/(0) + ^'(0)V e i 

3y 



(65) 



(66) 



and leave the precise choice of the value of m for later - it will 
be suggested by the asymptotic analysis. 
The corresponding inner equations are 



d,<& - -v„d p <& = 
e 

fy(V 2 <D,<S) = 



with 



MV ■ e m v4-<5MV 2< l>,<D), 

e z 

+2/(0) + ^-A'«D)V d , 
3y 



1 r / e 2 (VOV\»> 



(67) 



4/(<D) 



(V<D) 



(68) 



A. Leading order 

In the outer equations, g becomes the identity operator to 
leading order, i.e., O (O) (0 <O) ) = 1, and at the lowest order in e, 
we have 



V Z /'(0 (O) ) = 



(69) 



a Laplace equation that we know to be uniquely solvable for 
/'(0 (O) ) with Dirichlet boundary conditions at infinity. This 
boundary condition is even homogeneous (except possibly in 
a boundary region of order e), leading to the unique solution 
/'(0 (O) ) ee 0. This leaves us with the three possibilities ( °) = 
0, 1, of which (O) = or (O) = 1 are realized, according 
to the particular boundary condition on <O) . 



Again, = and 0=1 are solutions to the outer problem at 
all orders of e. Admittedly, the operator Q becomes indefinite 
at order e 2 for = and 0=1 [because of the denominator 
/(0)], but this does not matter, since the expression for dp 
alone is zero already at = and 0=1. 

The leading-order inner equation reads [g = 1 + (9(e)] 



1 - 



con 2 



4/(<DW>) 



dp(d pp ® m - 2/'(O (0) )) = . 



(70) 



Clearly, this is solved by O <0) = V 2 (1 - tanhp), which makes 
both the expression in brackets and in large parentheses van- 
ish. If we require m to be even, this solution is moreover 
unique (up to translations, which are eliminated by the re- 
quirement that the interface be at p = 0). For as soon as 
we assume (d),(, 0) ) 2 + 4/(<D (0) ), the mth power of the bracket 
expression will be positive, allowing us to use similar argu- 
ments as in Sec. IIII Cl between Eqs. J3 1 1 and J34l > to prove 
that dppO (0) - 2/'(<l> (0) ) = 0, and hence the bracket expres- 
sion must be zero, contrary to our assumption. Thus we do 
get a definite solution for <5 (0) from the inner equation, which 
moreover shows that at leading order of the inner expansion 
the second-order p derivatives of the operator V ■ Q'"V cancel 
each other. 



B. Next-to-leading order 

To simplify computations at the next order, we first expand 
V ■ Q m V up to formal order e°. This produces 



V ■ Q m V = 
1 

Vf 
1 



1 

•p 



1 1 „ e 2 (V(D) 2 \'" 



1 

vi' 



4/«D) 
1 



+ d p Vg^0> >a g^,p d p + O(e) . 



(71) 



Given that cD (0) is a function of p only, we realize that the third 
and fourth terms on the right hand side are (9(e), containing 
derivatives with respect to s or q of O, the fifth is even (9(e 2 ), 
so these terms may be dropped immediately in an expansion 
up to 0(1). The second term on the right hand side owes its 
existence the fact that Q is not exactly the projection operator 
on n [note that no such term is present in Eq. ( I59H and it has 
a prefactor of 1/e 2 due to the double derivative in p. This 
term which is desirable at leading order, because without it 
we would not have a determinate zeroth order solution <D (0) , 
is somewhat disturbing at the next order. Since the order of 
this term is 0(e"'~ 2 ), we can make it small by choosing m > 3, 
i.e., restricting ourselves to even m for the reasons discussed 
before, we set m = 4. Then the only remaining term on the 
right hand side of Eq. HH up to order e is the first term, 
which is the desired surface Laplacian. 
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Using this result, we can write the next-to-leading (nontriv- 
ial) order inner equation 

-v„3 p <£ (0) = M-^-d^g*%6fi m , 

S/i m = - jxo (1) + (*! + K2)dp& 0) 

_ ±hWV)Vf}, (72) 

again with X as given in Eq. ( 1421 . 

Note that we actually seem to have skipped orders here. The 
leading-order inner equation is formally (9(e~ 4 ), but once the 
zeroth-order inner solution is fixed, the differential operator 
V • g"'V is, according to dHJ, of order e max(0 > m - 2) only, so 
the order e 3 vanishes identically. The order e~ 2 is satisfied 
automatically, because the zeroth-order chemical potential is 
zero; the next nontrivial order is e -1 . Alternatively, one may 
say that the effective leading order has become e~ 2 . 

Keeping only order one terms of the surface Laplacian, we 
have V VI <5s> y/gg^dp = d ss + d qq , and the total linear operator 
in front of becomes —{d ss +d qq )£,. It is hermitian, because 
its hermitian factors commute. Hence, c? p O (0) is a left null 
eigenfunction. Multiplying (1721 from the left by it, integrat- 
ing with respect to p from -oo to oo, we obtain Eq. J50i . with 
M = M*y; to arrive at Eq. (|6), we have to set M = D s y/p s . 
Together with App. [C] this proves that the phase-field model 
of this section exhibiting a modified tensorial mobility has 
the correct asymptotic behaviour for small e, neither overcon- 
straining the system by adding, nor leaving it indeterminate 
by losing equations. 

VL CONCLUSIONS 

The intuitive approach to constructing a phase-field model 
for surface diffusion consists in using the chemical potential 
known from the nonconservative model to define a current, 
involving its gradient and a mobility that vanishes in the bulk 
phases, and to take the divergence of this current as the time 
derivative of the phase field. As has been shown in this arti- 
cle, this approach fails to produce the correct asymptotics in a 
subtle way. It does reproduce the equlibrium limit correctly. 

Next, the idea was explored that this failure might be 
remedied by making the mobility a tensor. After all, sur- 
face diffusion may be interpreted as highly anisotropic three- 
dimensional diffusion with a diffusion tensor that has zero 
eigenvalue in one direction. A straightforward realization of 
this idea fails in a rather drastic way, because restricting dif- 
fusion to the surfaces of constant phase field does not impose 
any functional dependence of <p in the normal direction. 

A modification of the tensorial mobility however leads to a 
model exhibiting the correct asymptotics. This model is given 
by the equations J64> . J66> with m — 4, 0, and (|9j- Larger 
even values of m also give valid models, a fact that might be 
useful in optimizing numerical performance. It is possible that 
m = 3 will lead to a valid model, too, but without positivity 
of the mobility tensor, it is difficult to prove that the inner 
equation ( 17 0> is solved by a unique phase-field profile. 



Since the condition (I52> destroying the asymptotic correct- 
ness of the scalar-mobility model is a strong restriction, it is 
conceivable that the model will satisfy the other equations 
while not being asymptotic yet and therefore keep some of 
its utility. It should however be emphasized that no positive 
assertions concerning its quantitative validity are available for 
this case. The problem of making reliable analytic statements 
about a nonasymptotic state is rather difficult. 

If the scalar-mobility model had a range of quantitative va- 
lidity, this would be restricted to some finite time interval, be- 
fore the full asymptotics kicks in. Even then, it may be argued 
that the modified tensorial model discussed here is probably 
more useful for numerical simulations. The point is the fol- 
lowing. In the scalar-mobility model, the inner equation de- 
termining the interface velocity appears only three orders in 
e after the leading order. A simulation must of course repre- 
sent the full model equations. Suppose we wish the interface 
velocity to be determined with an error not exceeding order e. 
Then the leading-order equation must be simulated with an ac- 
curacy <9(e 4 ) at least. If one takes a simple second-order accu- 
rate discretization for the gradient operators in the equations, 
the numerical error due to discretization alone will be 0(a 2 ) 
for a grid spacing a. To keep this smaller than e 4 , the grid 
spacing would have to scale with e 2 , which can be expected 
to lead to prohibitive computation times. Therefore, in the 
scalar mobility model, high-accuracy discretizations would be 
mandatory, even if the asymptotic error were controllable. 

On the other hand, in the model discussed in Sec.|V] as soon 
as the phase-field profile is represented with an error of order 
e or better, the effective leading order is only one order lower 
than the one determining the interface velocity, similar to the 
nonconservative case. Hence, reasonable accuracy should be 
attainable with grid spacings that scale as e, not as e 2 . None 
of the benefits of high-accuracy discretizations would be lost 
to the need of representing terms very accurately that are very 
small in the leading-order equation. 

Comparing the structure of the scalar-mobility model and 
the one based on a modified tensorial mobility, it is clear that 
the former model has much more intuitive appeal while the 
latter looks pretty complicated. But it has the virtue of pro- 
ducing the correct sharp-interface asymptotics. Very likely, in 
the future simpler and more elegant models will turn up doing 
the same thing. In the meantime, this one may be of use. 
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APPENDIX A: MATCHING CONDITIONS 

Let iff(x, y, z, t) - if/(r, s, q, f) be some arbitrary (sufficiently 
often differentiable) function of space and time obtained in 
solving the outer equations. We write the corresponding func- 
tion of the inner solution as ^(p, s, q, t), and suppress from 
now on, in this section, the dependence of functions on s, q, 
and t. Moreover, we write the coefficient functions in expan- 
sions with respect to e with simple subscripts indicating their 
order rather than superscripts in parentheses as in the main 
text. There, the notation is dictated by the fact that a subscript 
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would interfere with subscripts indicating tensor components; 
here, a superscript would interfere with the primes denoting 
derivatives. 

We must have the asymptotic relationship 

¥(p) ~ iff(r) = iff(ep) (p oo, e 0, ep 0) . (Al) 

Expanding both functions in powers of e, we get 



¥(p) = "i'o(p) + e^ 1 (p) + e 2 ^ 2 (p) + . 
t//(ep) = ^ Q (r) + etf/i(r) + e 2 ^ 2 (r) + .. 
= M0) + e[p^(0) + M0)] 



(A2) 



+ e z [p z -^(0) + p^(0) + ^ 2 (0)] + • ■ • , 

(A3) 

where the derivatives are to be taken for r — > +0, should they 
be discontinuous at r — 0. Analogous expressions with r — > 
-0 are obtained for the asymptotics as p — > -oo. 

Equating powers of e, we then successively get the asymp- 
totic relationships 



lim W (p) = M±0) , 

p— >±oo 

¥i(p) ~ p^ (±0) + ^(±0) (p^+oo), 



¥ 2 (p) ~ -pVo(±0)+P<A'i(±0) + ^2(±0) 



(p — > +oo) . 



(A4) 



(A5) 



(A6) 



Moreover, asymptotic relations such as ( IA5i can be decom- 
posed into statements about function limits 



lim d^iip) = ^(±0). 



lim o [>P 1 (p)-p^(±0)] = *Ai(+0) . 



(A7) 
(A8) 



APPENDIX B: LOCAL COORDINATES AND METRIC 

Once the coordinates s, q on the interface have been chosen 
and the third coordinate r has been specified via ( II It . the co- 
ordinate basis dl 2l i and, hence, the metric is determined. How- 
ever, in order to obtain explicit expressions and because it is 
desirable to interpret the metric in terms of familiar geometric 
concepts, we need a few additional considerations. 

For any curve r = R(s) on a surface, an accompanying or- 
thonormal frame can be defined consisting of the normal on 
the surface n at the considered curve point, its tangent t, pro- 
portional to d s R(s) (and equal to it, if s is the arclength), and 
the vector 1 = n x t. This triple, in the older German literature 
l27il sometimes referred to as the Bonnet-Kowalewski trihe- 
dron, is more commonly known under the name of Darboux 
frame 1 28] . It satisfies the following equations: 



(Bl) 



f d,n ) 




' K T\ 




n 


d s t 




—k a 




t 


, d s \ , 




K —T -a , 




{ 1 J 



where k is the normal curvature, i.e., the curvature of the pro- 
jection of the curve onto the plane spanned by its tangent and 
the surface normal, a is the geodesic curvature, i.e., the curva- 
ture of the projection of the curve onto the tangential plane of 
the surface, and r is the geodesic torsion. The total curvature 
k s of the space curve given by R{s) is related to the normal 
and geodesic curvatures via 



2 2 2 

k= a + k 



(B2) 



In the mathematical literature |27], the sign conventions for 
the normal curvature and geodesic torsion are opposite to the 
ones used here. 

If we identify our general curve on the surface with a coor- 
dinate curve of constant q from Sec. IIII Al then the Darboux 
frame will consist of the vectors n, tj, and t2 varying with s; if 
we identify it with a curve of constant s, the Darboux frame is 
given by n, t2, and — ti (it is defined right-handed). This pro- 
vides us with formulas for the derivatives of the normal and 
tangent vectors of the surface, 

dn 

ds 
<9ti 
ds 
dh 
ds 
dn 
dq 
d% 
dq 
dU 
dq 

from which we immediately get the expressions fl!2i for the 
vectors of the canonical coordinate basis {S n S s ,S q }. More- 
over, we obtain some simplifications. Because S a = d a r, we 
have the identities 



Kitl +Tlt 2 , 


(B3) 


-Kin + a\t2 , 


(B4) 


—Tin — aiti , 


(B5) 


K 2 tl - T 2 t] , 


(B6) 


-K 2 n - a 2 ti , 


(B7) 


T 2 n + a 2 t2 , 


(B8) 



de&a - d a &R 



(B9) 



These are automatically fulfilled, if a — r or /3 = r, but yield 
nontrivial relationships for a — s, (3 — q. 

Consider a point r = R(s,q) on the interface. Then S s = 
d s R(s, q) — ti, S q — dqR(s, q) = t 2 , because we have specified 
our coordinates to be arclengths (otherwise, £, and & q would 
not be equal, just proportional to the unit tangent vectors). 
From (IB5> and (IB8I . we deduce that 



r 2 n + ff2t2 = -Tin - citi 



(B10) 



Since the tangential vectors are linearly independent this 
implies a\ — a 2 — 0, i.e., the geodesic curvature vanishes and 
our coordinate curves are necessarily geodesies. Then their 
normal curvature becomes equal (up to a sign) to the spatial 
curvature, and formulas JB3I through JB8I reduce to the better 
known Frenet-Serret formulas I27l 12811. Secondly, we obtain 
from Eq. ( IB 101 that t\ = -t 2 , because the prefactors of the 
normal vector must match. 
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This means that the metric, Eq. H3\ . when restricted to the 
interface, will depend on just three parameters, viz. K\, k^, and 
t. But this is precisely the number of independent elements of 
a symmetric 2x2 tensor, so we can describe the most general 
metric this way locally. 

To complete the discussion of geometric relationships near 
the interface, let us note that formula (II 8i applied to A = & r 
at r - 0, yields 



K = V • II = K\ + K2 



(Bll) 



i.e., the mean curvature k is the sum of any two normal curva- 
tures corresponding to orthogonal directions, not just the sum 
of the two principal curvatures. This is the ultimate justifica- 
tion for using the same notation for the principal curvatures in 
Sec.|ll]and the normal curvatures in Sec. lIII Al 

Another interesting geometric quantity is the Gaussian cur- 
vature kq, given as the product of the two principal curvatures. 
A general formula in terms of the Darboux frame is 

k g = det (n, (ti • V)n, (t 2 • V)n) = n (d s ii x d q n\ . (B12) 

Evaluating this with the help of ( IB3I and ( IB6> . we find an 
expression for the Gaussian curvature in terms of normal cur- 
vatures and torsions 



KG = K\K2 - T . 



(B13) 



From this, one may conclude that the torsions have to vanish 
when k\ and K2 become the principal curvatures, i.e., take their 
extremal values. 



APPENDIX C: SHARP-INTERFACE LIMIT OF THE 
ELASTIC PART OF THE MODEL 

Loosely speaking, the title of this appendix refers to two 
things. On the one hand, it must be shown that Eq. (|9} re- 
produces the bulk elastic equations as well as the boundary 
conditions for the elastic problem. On the other hand, VT 
must be evaluated and shown to be the elastic contribution to 
the chemical potential of Eq. Q. 

The first task is easily accomplished. Inserting the outer 
solutions <p — 1 and = into (|9), we obtain 



j 1 
dp 



— . 



(CI) 
(C2) 



which are the mechanical equilibrium conditions in the solid 
and the liquid, respectively. The boundary conditions follow 
from the inner equation ( 1241 . which at leading order reduces 
to 



= d p (t$&) + <Ke) 



where we have used that the S y become the basis vectors of 
the Darboux frame at leading order and that these are inde- 
pendent of r and, hence, p. From ( IC3> . we have 



B E (0) - 8 S (0) - B E (0) - 



(C4) 



i.e., the three appearing components of the generalized stress 
tensor are independent of p. The matching conditions then 
provide us with 

<x„„ = lim 1<® = lim = -p , (C5) 

p— >— CO p— >co 

(C6) 
(C7) 



a nh = lim tfj = lim = , 

p— >— CO p— >co 

cr nt2 = lim W = lim W = , 



with the limit values in the nonsolid phase (p — > oo) follow- 
ing directly from (J9jl- At the next order, the phase-field equa- 
tions yield the capillary overpressure correction to \C5l . but 
we shall not consider this level of detail here. 

For the second task, we need the strain tensor 1A in curvi- 
linear coordinates to evaluate in the inner domain. To per- 
form the transformation, we first note that (S' a is the zth carte- 
sian component of S a ) 

1 / \ 8x' Bx' 

= 2^ + U ^W 
llBUiBx 1 BUj Bx> 
2\B^~da + da dp 



= 1 [dp(U0 a ) - UidpQ + BaiUj&j,) - UjB a & j p \ 

= \ [dpU a - Vdp6 a + daUp - Ud a 8p] , (C8) 



where some additional notations should be obvious: t/,- is a 
cartesian component of the displacement vector, U is the full 
vector. Introducing the connection coefficients or Christoffel 
symbols (of the second kind) 



(C9) 



= ndplW + tiflpfg? + t 2 3 p £< 1 ? + 0(e) , 



(C3) 



we can rewrite dC8> as 

= \ (Uafi + Up v ) , (CIO) 

where the semicolon in the last line denotes a covariant deriva- 
tive. This last result is evident - tensorial equations retain their 
invariant form when written with covariant derivatives. Since 
in a flat space the derivative with respect to a cartesian coor- 
dinate is automatically a covariant derivative, we could have 
written down the last equation from knowledge of the form of 
the small-strain tensor in cartesian coordinates. 

The fastest way to calculate the Christoffel symbols seems 
to be via their definition ( IC9I . There is an alternative formula 
using derivatives of the metric coefficients and known to most 
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physicists from general relativity courses, but while this pro- 
duces the same results, it does so via lengthy intermediate ex- 
pressions only. Since we will consider only the leading-order 
expression for V e i, we give the connection coefficients only 
to that order as well. At order e°, the metric Jl 3i becomes the 
unit tensor, so superficial consideration might suggest that, be- 
ing defined via derivatives of the metric, the Christoffel sym- 
bols should all vanish. However, one has to be careful to first 
take derivatives and then collect terms of order one, because 
derivatives with respect to r produce a factor 1/e when rewrit- 
ten in terms of the stretched coordinate p. Taking this into 
account, we find the following Christoffel symbols to generi- 
cally be nonzero at leading order: 



r 

ss 


= ~K\ , 




r 

11 


= -K2 , 




r 

sq 


= r = 

qs 


—T 


r s 

A rs 


= r s = 

sr 


Kl , 


r v 

rq 


1 qr 


T , 


r q 

1 rs 


_ p? _ 
— 1 sr 


T , 


1 rq 


= r« = 

A qr 


K2 



(Cll) 



All the remainig ones are either 0(e) or exactly zero. 
The components of the strain tensor are 

1 

rr U Yj p , 

U ss = U s , s + k x U r + 0(e) , 
Uqq = U m + K 2 U r + 0(e) , 



Uu r , s +^U s , p ] — A] L> - TU q + ()(£). 

Uur4 + -U 9 f]-K : U ll -TU, + 0(e). 



U sq = Uu s , q + Uq, s ) + TU r + 0(e). 



(C12) 



Since the left hand sides of these equations are of order 1, 
the terms having a factor 1/e can be of order 1 at most, too. 
This implies Jjfl = JjfJ, = uf^ = 0, hence the zeroth-order 
displacement components must not depend on p: 



Uf> = Uf\ S ,q). 
Uf> = Uf\s,q). 

Uf\s,q). 



uf 



(C13) 



From this and JC12L we can immediately conclude that sev- 
eral of the strain tensor components are independent of p at 
zeroth order as well: 



ii 

y sq 



UfXs,q): 
Uf q \s, q ). 
Uf q \s,q). 



(C14) 



The new feature in comparison with the two-dimensional case 
is that even a shear component of the strain tensor will be con- 
tinuous across the interface, i.e., in a simulation the nonsolid 



phase will possibly undergo a biaxial strain. This is not really 
a problem, keeping in mind that this phase is treated in simula- 
tions as if it were a solid phase with vanishing shear modulus. 
Shear strains in this phase do not correspond to anything phys- 
ical, and there is neither a shear stress nor energy associated 
with them. 

If we take limir^oo p = po as boundary condition for the 
pressure in the nonsolid bulk phase, then we know from 
Eq. (IC2> that the pressure will take that value everywhere, 
and we can integrate the equation for £ rr from dC4> with a 
known integration constant. Writing the result out in terms of 
the strains [using Q], we obtain, introducing the abbreviation 
h = ft(O (0) ): 



ho[2GUP+A(u«> + U«» + ug)] 
+ (l-h )~A(u«? + U«» + U q y)-p = 



-po 



(C15) 



We know the p dependence of all terms in this equation ex- 
cept uf^ , which suggests to solve for the latter quantity. This 
results in 



2Gh Q 



(2G + A- A)h + A) 



(C16) 



v 



The equations for £ ( vr and £' y , from dC4t imply together with 

the boundary conditions C^Jl and (£7} that \jf) = t/, (0) - n 
Then Eq. (|SJl turns into 



V 



(0) 



(C17) 



and we can conclude from dC 161 and dC14> that the elastic 
potential depends on p only via ho- Abbreviating t^vU^ = 
Ufj + Ufq and X = (2G + A - A)h + A, we get the compact 
form 



2Ghi\ 2Glii)A 



X 

\2 



X 2 



Vf =G(tr 2D U«») 

+ G [(u^f + (uf q y + 2(u^) 
According to ( ID9> and ( ID10L we have 



y(0) _ 



Jo 



(C18) 



(C19) 



so the only integral that actually necessitates a calculation 
when integrating dC18b is 



X 



1 (2Gh 2Gh A\ 



X X 2 
2G 

(2G + A- A) 2 j a 
2G 



Ja 



(X - A)(X + A) 

Y 2 



dX 



2G + A 



(C20) 
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all the other integrals reduce to integrating a constant from 
zero to one. 
We obtain 

v-~h D ^) 2 + ^(^) 2 . (C2i) 

afi 

Note that both terms are interface scalar invariants, taking the 
same form in any (2D) coordinate system. 

To relate this to the sharp-interface formula Q, we should 
transform it into an expression in terms of surface quantities, 
belonging to the outer equations. We have 

«2? = Uf?(s,q), 

-^iW^S- (02) 



,(()) 



the last equation following from dC16b by taking the limit p — > 

— oo. 

In the following, we drop the superscript (0). Hooke's law 
together with dC22> provides us with 



1 - V 

E 
1 + v 



(cr ss + cr qq - 2cr rr ) , 

(o- ss - (T qq j , (C23) 



and the prefactors are related to the Lame constants by (1 - 
v)/E = (2G + A)/[2G(2G + 3A)] and (1 + v)/E = 1/2G. 

Inserting dC23t into ( IC2U and renaming cr rr into cr m , it is 
straightforward to show that 

1 + v 



- (()) 1 -I- V / \ z 



>]'■ 



(C24) 



which is, up to a factor, of l/p s the elastic part of the chemical 
potential of Eq. Q, just written in curvilinear coordinates in- 
stead of cartesian ones. It is clear from this result that if we set 
[in Eq. j66H M = jM* = yD s /p s , Eq. d50b will become iden- 
tical to Eq. (|6} from the sharp-interface model, which com- 
pletes the proof that the phase-field model of Sec. [V] asymp- 
totically approaches the sharp-interface model of Sec.llTl 



APPENDIX D: USEFUL PROPERTIES OF THE PHASE 
FIELD FUNCTIONS 

In order to simplify it for the reader to find the actual re- 
lationships for the various functions involving the phase-field 
that are used in the text, they are collected here for reference 
(and concreteness). Often only certain properties but not the 
precise form of these functions is important. 

The chosen double-well potential is 



Its derivative is given by 

f'(<P) = 20(1 - 0)(1 - 74) . (D2) 
The "volume fraction function" h(<p) is 

h(cf>) = <p 2 (3 - 20) , (D3) 

having the derivative 

= 60(1 - 0) . (D4) 
Then a useful observation is that 



fUI» = I gA'(0) 



(D5) 



To solve the ordinary differential equation satisfied by the 
zeroth-order inner solution 



(D6) 



1 and 



dppO^ - 2/'((D (u, ) = , 

with boundary conditions lim p ^_o<, cD (0) (p) 
lirrip^co 4> (0) (p) = 0, we multiply by d p &°\ integrate and take 
the square root (with the correct sign) to obtain 

3 P <D (0) = -2O (0) (l - <D (0) ) = - l -h'(^) , (D7) 

which can be solved by separation of variables. The solution 
is, up to a translation in p, given by 



O (0) = ^(1 -tanhp). 



1 



(D8) 



Requiring the position of the interface to be at p — fixes the 
choice out of the one-parameter set of solutions, present due 
to the translational invariance of the differential equation. 

With the help of the second equality of (ID7> . it is easy to 
calculate certain integrals appearing in the asymptotic analy- 
sis. Those integrals typically contain the factor ((9 p <t> (0) ) ; to 
do the integral, it is then beneficial to replace one of the fac- 
tors (and only one) with -/!'(<E> (0) )/3. Integrals obtained this 
way have the structure 

r°° 2 

/= dpf(h(& 0) )) (<9 P 0> (0) ) 

U —CO 

= -- f dpf(h(¥ 0) ))h'(¥ Q) )d p ¥ 0) 

3 J —CO 

= - I J dO^ f{h(& 0) ))h'(& 0) ) = lf dh f(h 



(D9) 



This way, one arrives, for example, at 

r^fdp^ 1 -. (dio) 

xj —CO *J 



/(0) = 2 (l-0) 2 



(Dl) 
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